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Abstract 

Solvent-based thin-film deposition constitutes a popular class of fabrication strategies 
for manufacturing organic electronic devices like organic solar cells. All such solvent-based 
techniques usually involve preparing dilute blends of electron-donor and electron-acceptor 
materials dissolved in a volatile solvent. After some form of coating onto a substrate to 
form a thin film, the solvent evaporates. An initially homogeneous mixture separates into 
electron-acceptor rich and electron-donor rich regions as the solvent evaporates. Depending 
on the specifics of the blend, processing conditions, and substrate characteristics different 
morphologies are typically formed. Experimental evidence consistently confirms that the re- 
sultant morphology critically affects device performance. A computational framework that 
can predict morphology evolution can significantly augment experimental analysis. Such a 
framework will also allow high throughput analysis of the large phase space of processing pa- 
rameters, thus yielding considerable insight into the process-structure-property relationships 
governing organic solar cell behavior. 

In this paper, we formulate a computational framework to predict evolution of morphol- 
ogy during solvent-based fabrication of organic thin films. This is accomplished by developing 
a phase field-based model of evaporation-induced and substrate-induced phase- separation in 
ternary systems. This formulation allows most of the important physical phenomena af- 
fecting morphology evolution during fabrication to be naturally incorporated. We discuss 
the various numerical and computational challenges associated with a three dimensional, 
finite-element based, massively parallel implementation of this framework. This formula- 
tion allows, for the first time, to model three-dimensional nanomorphology evolution over 
large time spans on device scale domains. We illustrate this framework by investigating 
and quantifying the effect of various process and system variables on morphology evolution. 
We explore ways to control the morphology evolution by investigating different evaporation 
rates, blend ratios and interaction parameters between components. 

Keywords: phase separation, evaporation, Cahn-Hilliard equation, substrate patterning, 
organic solar cells, morphology evolution. 
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1. Introduction 



Organic solar cells (OSC) manufactured from organic blends [TJ [2] represent a promising 
low-cost, rapidly deployable strategy for harnessing solar energy [3j. In contrast to traditional 
silicon based solar cells, organic (or polymer) solar cells are low weight, printable on flexible 
substrate, and most importantly, can be produced at room temperature at very low cost. 

Solvent-based thin-film deposition technologies [4J (e.g., spin coating, drop casting, doc- 
tor blading, roll-to-roll manufacturing) are the most common organic photovoltaic man- 
ufacturing techniques. These techniques, especially doctor blading and roll manufactur- 
ing, are very attractive, due to their ease of scale-up for large commercial production. All 
solution-processing techniques usually involve preparing dilute solutions of electron-donor 
and electron- accepting materials in a volatile solvent. After some form of coating onto a 
substrate, the solvent evaporates. An initially homogeneous mixture separates into electron- 
accepting rich regions and electron-donor rich regions as the solvent evaporates (see Fig- 
ure [I]). Depending on the specifics of the polymer blend and processing conditions (spin 
coating time [5j E], solvent type [3 |8], nature of substrates), different morphologies are 
typically formed in the active layer. 
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Figure 1: Schematic illustration of phenomena during solvent-based fabrication techniques (left). Solution 
of donor, acceptor and solvent is coated on the chemically active substrate. As the solvent evaporates into 
the atmosphere, morphology evolution is triggered. (Right) An example thin film sandwiched between two 
electrodes. Internal structure represents the typical morphology of OSC, i.e. the mixture of electron-donor 
regions (black) and electron- acceptor regions (white) sandwiched between electrodes. 



The active layer of OSC is a blend of two types of materials: electron donating and 
electron accepting material. The active layer is sandwiched between electrodes (see Figure [j] 
right). The morphological distribution of the electron-donor and electron-acceptor subre- 
gions in the device strongly determines the power conversion efficiency (PCE) of OSCs 0[9]. 
In fact, every stage of the photovoltaic process is affected by this morphology. Consequently, 
there is immense interest to understand morphology evolution during fabrication [71 fTUl fTTj . 
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An important technological goal (that understanding morphology evolution will help achieve) 
is the ability to design fabrication processes to obtain tailored morphologies for high efficiency 
OSC devices. 

Current state-of-the-art approaches towards tailoring the manufacturing process are lim- 
ited to combinatorial trial- and-error-based experimental investigations UHl E2] • Further- 
more, the inability to experimentally visualize morphology evolution hinders the ability to 
quantify the effect of various process and system variables (such as evaporation rate, blend 
ratio, solvent type) on morphology evolution. Experimental techniques provide only limited 
data for analysis: (a) mainly limited to observations of the lateral organization of the top 
layer, and (b) mainly limited to the final morphology. In addition, visualizing 3D morphology 
remains challenging. These challenges serve as a compelling rationale for developing a com- 
putational framework that can model morphology evolution, thus significantly augmenting 
experimental analysis. 

Computational approaches to this problem exist, but are mostly limited to one scale: 
atomistic scale or macro scale. From a macro scale perspective the problem of thin film 
formation is well-studied. The series of work summarized in [13] link the macroscale film 
thickening during evaporation process with angular velocity of the coater, concentration, 
evaporation rate, and solution viscosity. However, morphology evolution is not analyzed in 
these studies. On the other end of the spectrum, morphology evolution in a typical OSC 
system was recently studied using molecular dynamics simulations [Hj. The authors were 
able to predict phase separation between two typical components in a cubical domain of size 
25 nm and only for 135 ns ) without incorporating the macroscale effects of evaporation or 
substrate. To the authors best knowledge, there exists no meso-scale approach that links 
morphology evolution at the nano scale with macro-scale phenomena like evaporation, and 
substrate patterning. The development of such a model will provide significant advantages, 
it will in particular: 

• Serve as a powerful tool to analyze morphology evolution over time in three dimensions. 
This can be used as a "stereological microscope" to visualize morphology evolution from 
early stages until the formation of the stable morphology. It is worth mentioning that 
three dimensional experimental reconstruction of the final morphology is possible using 
electron tomography [9J [15l HE] but this approach to polymeric systems is exceedingly 
rare because of the required proper contrast between components. 

• Allow for independent control over various process and system variables, thus making 
it easy to isolate factors affecting the process (such as substrate patterning, solvent 
annealing, or blend ratio). 

• Ability to perform high throughput analysis. Such an ability allows automated ex- 
ploration of the phase space of various manufacturing and system variables (such as 
different blend ratios, evaporation rates, solvent choices, effect of substrates) to un- 
derstand their effect on morphology and performance. This opens up the possibility 
of data-driven knowledge discovery to understand the effects of different competing 
phenomena and, subsequently, tailoring the manufacturing process. 
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In this work, we develop a computational framework to model morphology evolution 
during fabrication of organic solar cells. We formulate a model that takes into account all the 
important processes occurring during solvent-based fabrication of OSCs. The model is based 
on a phase field approach to describe the behavior of multicomponent system with various 
driving forces. We develop an efficient numerical framework that enables three dimensional, 
long time-scale analysis of the fabrication. We illustrate the framework by investigating the 
effect of independent control over various external conditions: solvent evaporation, blend 
ratio, and interaction parameters. We further quantify the interplay between the solvent 
evaporation and diffusion within the film. To the best knowledge of the authors, this is the 
first comprehensive effort to construct a virtual framework to study 3D morphology evolution 
during solvent-based fabrication of organic solar cells. 

2. Physical phenomena affecting morphology evolution during fabrication 

All solution-processing techniques usually involve preparing dilute solutions of electron- 
donor and electron-accepting materials in a volatile solvent. After some form of coating 
onto a substrate, the solvent evaporates. An initially homogeneous mixture separates into 
electron-accepting rich regions and electron-donor rich regions as the solvent evaporates. 
Depending on the specifics of the blend and processing conditions (spin coating time [5j [6], 
annealing time [TTJ, HE] , solvent type [3 [8] , nature of substrates) , different morphologies are 
typically formed. The two materials usually used in fabricating OSC are a conjugated poly- 
mer and a fullerene derivative, The conjugated polymer is the electron donor material, while 
the fullerene derivative is the electron-acceptor. We will use the terms electron- accepting 
material (electron-donating material), acceptor (donor) and fullerene (polymer) interchange- 
ably in this article. 

There is a rich and complex collection of interacting phenomena that direct the mor- 
phology evolution during solvent-based fabrication. Phase-separation [6J [HI EUJ [21] is a key 
phenomena triggered by the evaporation of the volatile solvent. The atmosphere on the free 
surface determines the evaporation rate of the solvent. The resulting morphology can form 
multiple phases, making the system a multi- component and multi-phase system. In addition, 
the morphology on the substrate's surface may differ from that in the 'bulk' state j5j [22] . 
Both free surface and substrate influence the organization of the morphology and directly 
affect characteristics of the devices. In particular, chemical and physical patterning of the 
substrate have been shown to affect morphology evolution. In order to enable predictive 
modeling, each of the following phenomena must be included in the computational model. 

2.1. Evaporation induced morphology evolution 

The rate at which solvent is removed from the top layer depends on various factors 
(e.g. solvent volatility, spinning velocity during spin-coating). During evaporation, the 
initially dilute solution becomes enriched in the two solutes due to depletion of solvent. This 
enrichment results in increased interaction between the solutes and triggers morphology 
evolution. The evolution is critically determined by the evaporation profile and diffusion of 
solutes (and solvent) within the film. 
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Figure 2: Scheme of typical ternary phase diagram. (Left) Typical characteristic curves: binodal line (solid 
line) and spinodal line (dotted line). The spinodal line is related to the position where: q^q^z — d ^ = 0. 
The binodal curve is related to the equilibrium phase boundary between the single phase and the phase 
separated region. During typical fabrication, initial dilute solution, (jP s , is pushed into two-phase region 
(delimited by the spinodal line) as solvent evaporates (cjP s — > <j>\). Solution separated into polymer-rich ((//) 
and fullerene-rich phase (</>")• The equilibrium compositions of two phases lie on the binodal curve. Line 
connecting corresponding equilibrium compositions are called tie-lines (right). 

As the solvent is removed, the blend is pushed into the spinodal range (immiscible con- 
ditions) and induces phase separation (see Figure [2] left). Phase separation (or spinodal 
decomposition) is a mechanism by which solution separates to create phases of different 
properties. Under these conditions, the solution is unstable and even small fluctuations lead 
to fast phase separation. The composition of phases changes and reaching equilibrium con- 
centrations. Usually, the solution separates into two phases. One phase is rich in donor, the 
other phase is rich in acceptor material. The exact composition of phases is determined by 
thermodynamic conditions [23]. The creation of the two phases is followed by slow coarsen- 
ing. The kinetics of phase separation and coarsening is affected by the kinetics of solvent 
removal. 

2.2. Substrate induced morphology evolution 

In a confined thin film geometry, the evolution of the morphology close to the walls/interfaces 
can be significantly different from that of the bulk. In particular, chemical interactions be- 
tween the substrate and the solute as well as surface patterning can significantly affect 
morphology evolution. Recent experimental studies of OSC explore this possibility [5j [24] . 
By changing substrate properties and by (chemically or physically) patterning the substrate, 
it was possible to direct vertical segregation, percolation and control phase separation. 

In certain systems, crystallization is an additional mechanism of morphology evolution. 
While incorporating crystallization is relatively straightforward in the current framework, 
we primarily focus on modeling evaporation and substrate induced spinodal decomposition 
in this paper. This is in line with experimental results [25] which seem to suggest that phase 
separation is a key mechanism in morphology evolution for polymer-polymer blends. 
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3. Phase field approach to model mult i- component evolution 

Phase field methods have been used to model morphology evolution in heterogeneous ma- 
terials typically consisting of grains or domains characterized by different structure, orienta- 
tion or chemical composition. These methods are highly versatile (due to a diffuse-interface 
formulation) and easily represent the evolution of complex morphologies with no assumption 
made about shape or distribution of domains. Various thermodynamic driving forces for 
morphology evolution, e.g. bulk energy of the system, interfacial energy, substrate energy 
can be easily introduced. Additionally, the effect of different transport processes, such as 
mass diffusion, convection, or shear rates, can be directly introduced and exploited using 
this technique. 

The advantages and flexibility offered by the phase field method provide an ideal frame- 
work to model morphology evolution during fabrication of OSC. In particular, the mecha- 
nisms that direct morphology evolution (spinodal decomposition and crystallization) can be 
naturally modeled using this method. They have been well studied for alloy system (den- 
dritic growth and spinodal decomposition in alloy system [261 EZj)- Additionally, various 
driving forces for phase transformation, such as effect of substrate and evaporation, can be 
introduced in a very straightforward way without substantial model reformulation. Finally, 
phase field methods can be scaled to predict morphology evolution for device-scale problems, 
which is currently impossible using any other framework, like molecular dynamics. 

4. Mathematical model of evaporation— and substrate— induced phase separation 

In this section, we formulate the phase field model to simulate morphology evolution in a 
ternary system with solvent evaporation and substrate interaction included into the model. 

4-1- Ternary phase field model 

Formulating a phase field model usually consists of three stages: (1) identifying order 
parameters; (2) postulating free energy functional that depends on the order parameters; 
(3) constructing governing equations that describe the evolution of the system towards a 
minimum energy state. 

4-1.1. Phase field order parameters 

We consider a ternary system consisting of polymer, fullerene, and solvent. We denote 
the volume fractions of polymer, fullerene and solvent as (j)f and S , respectively. We 
set the volume fractions & as the conserved order variables (since <j) p + (j)f + (j) s = 1). Note 
that during evaporation, while volume is not conserved, the sum of the volume fractions - 
by definition - is conserved. 

4.1.2. Free energy functional 

In the second step, we construct the energy functional for this system. Energy, F, consists 
of two bulk terms: homogeneous energy, and interfacial energy between phases. The total 
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energy is given as: 



F((/> p ,(f>f,(f> s ) = 



v 



=p,f 



dV + F s (4> p ,4> f ,x) 



(1) 



The energy of the homogeneous system, /, also called configurational energy or free energy 
of mixing, is the quantity which governs phase separation. Homogeneous energy depends 
only on local volume fractions and is at least double- welled. Wells correspond to the equilib- 
rium concentrations of separated phases. In contrast, the interfacial energy depends on the 
composition gradient and is scaled by an interfacial parameter e (see second term in Eq. [j]) . 



In the current work, we assum^|e p = 6f = e. 

We construct homogeneous energy, /, using the Flory-Huggins formulation [23l l30l I3T] . 
which is suitable for polymer solutions. ^\ According to this theory, the energy of the system 
is given by: 



f((f>p,(f>f,(f>8) = 



RT 

V 7 



N, 



(2) 

where: R is the gas constant; T is the temperature; <^ is the volume fraction of component 
z; Ni{= V 1 /V r ) is the degree of polymerization of component i. Here, V 1 is the molar 
volume; V r is a reference molar volume (e.g. solvent); and Xij is the Flory-Huggins binary 
interaction parameter between component i and j. The interaction parameters and degree 
of polymerization define the shape of the binodal and spinodal lines (Figure [2]). 

One of the advantages of the phase-field method is the ability to introduce additional 
driving forces to the model in a straightforward way. For instance, the effect of the substrate 
is a surface energy term, F s ((j) p) 0/, x) that is simply added to the total energy of the system 
(Eq.[T]). A generalized form of the substrate energy is given in Eq.|3j Patterning is introduced 
through the space dependent functions p p (x) and p/(x) in Eq. [3j These functions determine 
the geometry of patterning. At a point x on the substrate that is chemically tuned to 
component z, the value of pi = 1, otherwise pi = 0. The chemical specificity of patterning 
i is reflected in the parameters ji l and h\ Parameter \i l is a chemical potential favoring 
component i at the substrate, and h p expresses the way interactions between the components 
near the substrate are modified by the presence of a pattern at the substrate j30j [32] . 



y ps ' 



RT 

V 7 



I 



f s dS 



RT 

'V 7 



Js 



+ 



D+PfWfaftf + hffi)] ds 



(3) 



1 Interfaces between different phases of different composition have different interfacial energies. However 
due to limited knowledge about interface properties, this simplification is common [28) 129]. 

2 Following [29 , we include an additional term, to the Flory-Huggins energy function. This 

modification allows to improve the efficiency of computation, while not changing the rate of morphology 
evolution or the shape of forming structure. The parameter j3 is set to a small value 0.001 RT /V s . 
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4.1.3. Governing equations 

Once the energy of the system is specified, the governing equations of the evolution can 
be formulated. This is usually done by defining the chemical potentials of the system. The 
chemical potential quantifies how much the energy changes when the configuration changes. 
The chemical potential for polymer and fullerene are \± v = 8F/8(j) p and jif = 5F/5cj)f, 
respectively. Next, using Fick's First Law for the flux (J p = — M p V ii p ) and the continuity 
equation (d(j) p /dt+ V J p = 0) we get the governing equation for each component. We consider 
only two of the three variables as independent (since (j) s = 1 — (j) p — (j)f). The resulting Cahn- 
Hilliard-Cook equations are given by: 



d(j) p h d(j) p 

— — -\- ii — 

dt h curr dh 



Jp v \ <% 



M„V [ --Z":^' - e 2 VV P + (1 - H (h)) 



df s {4>p,(f>f) 
d(j> p 



+ 4 (4) 



d4> f h d(j> f 
dt h curr dh 



\ 9(f) f d(f)f 



<S (5) 



where Mi = D / \d 2 j \deai I d 2 4>i) is the mobility of component i. The diffusivity, D, is a linear 
combination of self-diffusivities of all components and their volume fractions: D = D p (j) p + 
Df(f)f + D s (j) s . The energy of ideal solution, f ideal, 18 used to link mobility with diffusivity and 
to comply with classic Fick's law. The ideal solution is one with zero interaction parameters. 

The advection term (second term in LHS of both equations) accounts for the change in 
height of the film and is scaled according to the height. Scaling term h/h curr is zero at the 
bottom surface (h = 0) and one at the top surface (h = h curr ) ) where h curr is the total current 
height of the film. Substrate effects are included in the third RHS term of both equations. 
This term enters the equation only for the bottom surface (zero height h — 0) where the 
Heaviside function, H{h) = 0. The last term of RHS (in both equations) is the Langevin force 
term, t; p and This term mimics the conserved noise due to fluctuations in composition. 
We set the noise to be a Gaussian space-time white noise with the following constrains: 
xi)) = and (&(*i,xi)&(* 2j x 2 )) = 2M l RT/V s 5{t 1 - t 2 )V 2 5( Xl - x 2 ), where 5 is the 
Kronecker delta. Variance is determined by the fluctuation-dissipation theorem (FDT) [33] . 
Cahn-Hilliard equation with noise considered is called the Cahn-Hilliard-Cook equation [3^] . 

This formulation allows for a natural extension to include crystallization. To do this, an 
additional phase variable must be considered, the energy functional expanded accordingly 
and the governing equation formulated. 

4-2. Modeling solvent evaporation 

From a mathematical perspective, an evaporation process is classified as a moving bound- 
ary problem, since the height of the thin film changes over time (h = h(t)). Following [35] . 
we explicitly trace the surface evolution. We assume that solvent is the only component that 
evaporates from the top surface. Solvent lost from the top layer results in height decrease. 
The rate at which the height decreases, u = dh/dt, is proportional to the flux of the solvent 
out of system, 
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where J® is the normal component of molar flux (of solvent into the air). This equation 
constitutes the mass balance for moving film-air interface: u((j)^ — </>f) = V s (J® — Jf) with 
two assumptions: only solvent evaporates ( = J J = 0) and solvent content in the air is 
very low (cj) a s <C 1). The superscript a and s denote air and film, respectively. The molar 
flux of solvent can be further linked with evaporation rate of the solvent, k ei and the average 
content of the solvent at the top layer, 0^ op , using equation: 

v s r s = k e ft s op . (7) 

In this work, we assume that the solvent evaporates uniformly from the top surface and film 
height decreases homogenously. This assumption simplifies analysis and quantification of 
the competing effects of evaporation and substrate shown in the results section. We investi- 
gate the effect of inhomogeneous evaporation on the film surface evolution in a forthcoming 
publication. 

The evaporation rate, fc e , depends on various parameters including the solvent par- 
tial pressure, solvent vapor pressure, temperature, air flow and is specific to the fabri- 
cation process. For spin-coating, the evaporation rate depends on the angular velocity 
k e = k^uj 1 / 2 ) [36J, while for solvent annealing or drop-casting, the evaporation rate can 
be assumed constant (k e = const). We apply boundary conditions at the top surface to sat- 
isfy the balance of two other components within the film. Removal of the solvent from the 
top layer, results in enrichment of the volume fraction of polymer and fullerene. Therefore, 
to account for this enrichment we apply Neumann boundary conditions at top surface for 
two other components: J^ s = —J^cj)^. ^ 

4-3. Challenges 

Although the phase field method is a well used technique for simulating the morpho- 
logical evolution of a wide variety of materials and processes [371 EH EH] , employing it to 
predict morphology evolution during fabrication of active layer for organic solar cells requires 
resolving several challenges. We detail these challenges below. 

4.3.1. Physical Modeling Challenges 

Morphology evolution during fabrication of organic solar cells is an intrinsically multi- 
scale process both in time and space which makes this process difficult to solve accurately 
and efficiently using reasonable computational resources. The governing equations (Eqns. [IJ 
[5]) are forth-order nonlinear partial differential equations, which are difficult to solve numer- 
ically. The complexity of these equations is related to two competitive subprocesses: phase 
separation and coarsening. These processes occur at two widely different temporal and spa- 
tial scales, all of which must be resolved properly. Initially, very fast phase separation is the 
dominant process and followed by slow coarsening. Phase separation is a fast process and 



3 Applied boundary conditions guarantee that content of polymer and fullerene is conserved, that was 
confirmed by thorough tests performed. 
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results in thin layer creation. In contrast, coarsening is slow process consisting of rare events 
and involving merging of bulky regions. 

Another numerical challenge is related to the evaporation process. In case of the solvent- 
based thin-film deposition, the volume fraction of solvent changes from 99% to almost zero. 
This poses a severe problem for numerical techniques that have to reliably model the huge 
change in domain size in 3D. 

4-3.2. Computational Challenges 

A key objective of the formulation is the ability to predict morphology evolution at the 
device scale. The inherent complexity of the process makes this objective demanding. This is 
because we are interested in resolving nano scale morphological evolution while investigating 
device-scale domains. From a computational perspective this involves solving differential 
equations with a very large number of degrees of freedom (3> 10 7 ), which cannot be solved 
using current serial processing machines. This require developing modules heavily based on 
parallel processing, including applying domain decomposition strategies [391 140] . 

4.3.3. Materials Science Challenges 

The phase field method is a generic technique, and can be applied to almost any type of 
system undergoing morphology evolution. Thus, in order to provide quantitative prediction, 
it is necessary to determine material-specific set of parameters, both thermodynamic and 
kinetic. These parameters very often show compositional, directional or temperature depen- 
dence, which pose additional difficulties. In this regard, significant work has been done to 
determine several parameters for materials and systems utilized for photovoltaic applications 
using both molecular dynamics j4TJ [14] and experimental techniques [42l [25l [43] . 



5. Numerical model 

5.1. Strong form of the split Cahn-Hilliard equations with boundary conditions 

The strong form is formulated as follows: find (j) p) (j)f : fix [0,T] — >► K. and \i v ^\if : 
fix [0, T] —> R (/x is an auxiliary variable) such that: 



h 

u- 



dt h curr dh 
d</>f h <% 

dt h curr dh 



lh = jnr- £ v % + (1 " H(h)) 
H f = ^- e 2 V 2 0/ + (1 - H(h)) &fs 



d<j)f d<j)f 

V(M p fip) ■n = h IXp on 

\7(Mffj,f)-n = hp f on 

p (x, 0) = <£j(x) in 

^(x,O) = 0j(x) m 



n x [o,t], 


(8) 


x [0,T], 


(9) 


ft x [0,T], 


(10) 


ft x [0,T], 


(11) 


r t x[o,T], 


(12) 


r t x[o,T], 


(13) 


ft, 


(14) 
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(15) 
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We apply Neumann boundary conditions on the top surface, to account for polymer and 
fullerene fraction increase at the top layer due to solvent evaporation. Flux of polymer 
equals to — J^0 P , while flux of the fullerene h jlf equals to —J^cpf. On other boundaries, we 
apply zero flux conditions for both variables: volume fraction, and chemical potential, 

Mi- 

5.2. Domain rescaling due to evaporation 

To account for domain change due to solvent removal from the top layer, we introduce a 
coordinate transformation j44j [45] : 

# = (16) 

where h is the height coordinate scaled according to current total height of the film, h curr . 
This coordinate transformation permits recasting the model equations into a system of equa- 
tions with fixed boundaries which is more convenient for numerical solution. In this way, 
there is no need for remeshing. 

Since we assume homogeneous evaporation, film has an uniform height h curr {t) ) at any 



point in time, t. Recasting the model equation using the coordinate transformation [16 
involves defining the gradient operator in the new coordinate system: 

V = e ^m +ei te i (17) 

where the first component is the direction along the height of the film. Strong forms of the 
recast equations are given by: 



d(f) p t $ d(j), 
~dt 

d(j)f d d(f)f 

~dt 1 " 



+ U h^lf& = V-(M p V//p)+& infix [0,7], (18) 

+ U J^W = ^-( M /^/)+0 into*[0,n (19) 

Mp = ^-^V% + (1-H(h))?£- mfix[0,T], (20) 

H = §f-^% + (l-H(h))^ mfix[0,T], (21) 

V(M^)-n = onT t x[0,T], (22) 

V(M f fi f ) • n = on T t x [0,T], (23) 



The advection term (second LHS term of Eqns 18 and [19]) is non zero only along the height 
direction. Velocity corresponding to this term, u : is proportional to the molar flux of the 
solvent out of system: u = — V s J* = k e ^l° p ( Eq. [fj]). 

5.3. Spatial discretization of governing equation 



We use the finite element method to solve the governing equations (Eqns [18]- 23). We 
solve equations in the split form, to avoid constraints related to the continuity of the basis 
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functions when using the finite element method with a primal variational formulation [46] . 
More precisely, the standard C°-continuous finite element formulation is not suitable for 
forth-order operators, and consequently basis functions which are piecewise smooth and C 1 - 
continuous should be utilized |47l [48] . However, there are only a limited number of finite 
elements that fulfill the above continuity condition, especially in two and three dimensions. 
The weak form of the split Cahn-Hilliard equation is given by: 

O, (/> Pjt ) Q + O, Ui$/h curr (fi p ^ n + a (w, M p p p ) Q + a (w, p p ) n = (w, h^) 

K <f>f,t) n + K ^/h curr 0u)n + a K M fVf)n + a K 
df \ ( 2 v ( dfs 

9/ A , N , , 2jL N , ( df a 





= (w. 




(24) 




= (w 




(25) 


> y 


I- 





(26) 



O, /i/)^ + a (w, e 2 ; )^ + ^, |^(&)) = ( 27 ) 



where n; G are weighting functions, (•, •) is the L 2 inner product on fi, a(-, •) is the 

energy inner product on /i^. and define natural boundary conditions. Second term in 
LHS of Eqns [24] and [25] accounts for domain size change due to evaporation. Fourth term 



in LHS of Eqns [24] and 25 accounts for the conserved noise, where pi is the vector of the 
stochastic flux terms. Fourth term in LHS of Eqns [26] and [27] accounts for substrate effect 
and are included only for surface elements belonging to bottom boundary T^. We note that 
substrate term is not a typical boundary condition but an additional term resulting from 
the additional energy in the system. 

We use the Galerkin approximation to solve the two split Cahn-Hilliard equations. We 
define (j) p G S p and G S p to be the finite dimensional approximation of polymer and 
fullerene volume fraction fields, fi p G A4 p jJj G M 1 } to be the finite dimensional approx- 
imation of polymer and fullerene chemical potential fields and w h G V h to be the finite 
dimensional weighting function. The approximate solutions are computed on the following 
function spaces: 

S h p = {4> h p \4> h p e H 1 ^), 4> h p e P k (W) Ve} (28) 
K = {/i>J G H\n), 4 e P k (n e ) Ve} (29) 

S h f = {(f> h f \(t> h f e H\D), 4> h f e P k (tt e ) Ve} (30) 
M) = e H^Q), \x) e P k (n e ) Ve} (31) 

V h = {w h \w h e w h e P k (Q e ) Ve} (32) 

with P k (Q e ) being the space of the standard polynomial finite element shape functions 
on element f2 e , where k is the polynomial order. We additionally introduce the SUPG 



stabilization term for the advection terms in Eqs 24 and 25 |49| 
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We discretize the conserved noise by generating Gaussian- distributed random numbers 
for each component of flux p that satisfy (pi(t\,x\)) = and (fopy) = 1. We generate Tri- 
dimensional vector (n s d = 2,3 depending on 2D or 3D problem) for each component (polymer 
and fullerene) in each node at each time step, p pi = g p pi. Parameters g p = ^J2M P RT / (V s At) 
and q = ^/2M f RT/(V s At) (where At is size of the time step) account for fluctuation- 
dissipation theorem requirement. For more details see [501 [27] . 

We use linear basis functions for variables (j) p) (j)f ) p p) pf, p p and pf. We also tested 
quadratic and cubic basis functions but have not noticed any significant improvements. This 
choice is based on an extensive analysis of the Cahn-Hilliard equation in [5T] . 

5.4- Temporal discretization of governing equation 

In our approach, we take advantage of implicit time schemes due to their unconditional 
stability. Explicit methods are often intractable due to severe restrictions on size of time 
step (~ Ax 4 ) arising from the stiffness of the equation. This makes them computationally 
prohibitive even for simple problems. Consequently, implicit methods arise as a natural 
alternative. Since they are unconditionally stable they allow for much larger time step. 
However, because of nonlinear nature of the Cahn-Hilliard equation, implicit schemes require 
nonlinear solvers. The popular implicit time schemes in this context are Euler Backward 
and Crank-Nicholson methods [52l l53l 147] . 

Because of the multiple temporal scales that occur during phase separation, we tested 
adaptive various time stepping strategies [5T] . We noticed significant improvement in terms 
of number of times steps required to reach steady state as well as in total run times [5T] . 
In these time stepping strategies, step size is adjusted on the basis of the error between 
solutions of different order. However, whenever noise is considered, the error computed 
using such strategies is highly affected by the noise and cannot be used in standard time 
stepping strategies. Therefore, we use a Euler Backward scheme with a heuristic strategy to 



adjust the time step as used in [53] (see Appendix A). 



6. Technical details 

To solve the nonlinear system of two split Cahn-Hilliard equations, the formulation is 
linearized consistently and a Newton-Raphson scheme is used. To solve large problems with 
several millions of degrees of freedom, we use a domain-decomposition based mesh-partitioner 
to divide the mesh and distribute it across computational nodes. In our framework we use 
the ParMETIS partitioner [39j HQ] . Additionally, to enable parallel solution of the algebraic 
systems, we use the PETSc solver library [5H [55]. All results reported are obtained using 
the Generalized Minimal Residual Method. 



7. Results 

In this section, we showcase the formulation by investigating the effect of evaporation 
rate, blend ratio, degree of polymerization as well as effect of choice of solvent on morphology 
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evolution. Finally, we investigate the possibility of additional control of morphology through 
substrate patterning. 

We investigate a representative OSC system. Such a system consists of polymer, fullerene 
and solvent. In the default configuration, we consider degree of polymerization of solvent 
N p = 5, fullerene Nf = 5 and solvent N s — 1. Degree of polymerization is computed 
on the basis of molar volumes of the components. We assume V s = 100 cm 3 /mol, V? = 
500 cm 3 /mol, V p = 500 cm 3 /mol. Interaction parameter between polymer and fullerene 
reflects the low solubility of two components and is set asx s / = 1.0. Interaction parameters 
between solvent and fullerene or polymer are assumed to be much lower: Xfs — 0.3, Xps = 



0.3. In subsection |7.3[ we investigate the effect of varying degree of polymerization to 
N p = 20 and N p = 100. In subsection 7.4, we investigate the effect of solvent and various 
interaction parameters. All these values are representative for OSC and correspond well 
with experimentally determined values [25J. We take the solvent self-diffusion coefficient 
as D s = 10 _9 m 2 /s, which are typical for solvents (chlorobenzene, chloroform and xylene) 
used in OSC [56]. Self-diffusion of polymer and fullerene is much lower. We assume that 
D p = D f = 10~ 3 D S . 

We present one, two and three dimensional results. For each simulation we generate one 
mesh in the reference coordinate system. Mesh density is adjusted to the estimated width 
of the interface between polymer and fullerene. We use the interface width as a metric 
to determine mesh density to accurately capture the dynamics of the interface evolution. 
Our detailed analysis presented in [5]] showed that at least four elements per interface are 
required to capture the physics of phase separation and coarsening accurately. 

In [51], we showed that the analytical estimation of interface width is fairly accurate for 
two and three dimensional cases, even though it was derived for the one dimensional case 
in [26] . The interface width is defined as a distance required to span by the concentration pro- 
file across the accessible range of the composition: 5 = A(f) e /(d(f)/dx)\ ( f )c = A0 e ^/e 2 /A/ max . 
For more details see [26], EE]. The profile of concentration, and subsequently width of the 
interface, depends on the interfacial parameter e, interaction parameters, degree of polymer- 
ization. In ternary system range of the composition changes with time at the intermediate 
stages. Consequently, the interface width also changes with evaporation. To estimate inter- 
face width, we consider the case when the interface width is the smallest. This corresponds 
to the fully evaporated, binary fullerene-polymer system. For the system modeled in this 
paper, the interface width between polymer and fullerene varies from 12 nm (for N p = 5) to 
Gnm (for N p = 100). 

The interfacial parameter e is also estimated for a binary system consisting of polymer 
and fullerene. Following the analysis in [26, [29], one can link e with the interfacial energy, 
7 P /, between polymer-rich and fullerene-rich phases. We assume j p f = 33 mJ/m 2 , which 
is comparable with interfacial energies for organic compounds. For more details regarding 
detailed analysis and computations, we refer the reader to [29]. Computed value of e 2 is 
3.57-10- 10 J/ra (for N p = 5), 2.62- lO" 10 J/m (N p = 20) and e 2 = 2.05- lO" 10 J/m (N p = 100). 

The initial solvent fraction is chosen to guarantee that the ternary solution is homoge- 
neous and there is no phase separation at time t = 0. In 3D simulations, we start with 
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a solvent fraction (f) s = 0.66. In 2D simulations, we start with a solvent fraction 0.75. In 
ID simulations, we start with a solvent fraction 0.6. The higher initial volume fraction of 
solvent in 2D simulations is required to cover the wide range of degree of polymerization 
and interaction parameters investigated. Simulation is stopped when the fraction of solvent 
within the film is 0.05. At this time, the diffusion coefficient reduces significantly, and the 
morphology is frozen. 

In two dimensional simulation, we generate meshes that consist of 250 x 100 linear ele- 
ments. The computational domain is rectangular with height L y = 1, and width L x = 2.5. 
Height changes from Ly = 1 to 0.26, when all the solvent evaporates. In three dimensional 
simulations, we generated meshes that consist of 230 x 230 x 70 linear elements to discretize 
a computational domain that is a rectangular prism of length L x = 3.3, breadth L y = 3.3 
and height L z = 1.0. Height changes from L z = 1 to 0.34, as the solvent evaporates. In one 
dimensional simulations, we generate meshes that consist of 100 linear elements to discretize 
a height L x = 1. Height changes from 1.0 to 0.42. 

The large number of degree of freedom (~15 million for 3D simulation) require using 
parallel solvers and domain decomposition. The average run time for a three dimensional 
case is around 50h solved using 256CPUs. The average run time for two dimensional case 
is around 1.5h using 8CPUs. Detailed scalability analysis of the solver has been reported 
elsewhere I5T1. 



7.1. Influence of evaporation rate 

Evaporation of the solvent from the top surface is one of the key phenomena that induces 
phase separation. The way morphology evolves is an interplay between two competing 
processes: (i) evaporation of the solvent from the top surface, and (ii) the diffusion of the 
solvent within the film. This interplay can be expressed as a mass Biot number. Biot number 
is defined in Eq. [33] and expresses the ratio between external mass transport by evaporation 
of the solvent from the top surface and internal mass transfer to the top surface by diffusion. 

* - ^ (33) 

where L is the characteristic length - height of the film, D is the diffusion coefficient of the 
thin film and k e is the mass transfer coefficient - evaporation rate. We compute the Biot 
number for initial height and solvent diffusion coefficient. In such case, the solution consists 
mostly of solvent, and the diffusion is dominated by it. 

To better understand the interplay between evaporation and diffusion, we first perform 
experiments in ID. We consider the default case of system variables ( N p = Nf — 5, N s — 1, 
X P f = 1, Xps = Xfs = 0.3). In Figure |3j we plot change of height for three Biot numbers as 
a function of time. In the same figure we also show the ID volume fraction profiles of one 
component, polymer, as function of height. Initially, polymer is distributed uniformly along 
the height and its volume fraction is 0.2. Similarly, the initial volume fraction of fullerene 
is 0.2, since the blend ratio between polymer and fullerene 1:1. For the symmetric case, 
the fullerene profiles mirror images of polymer profiles (and hence are not plotted). We 
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investigate the cases with three Biot number: equal to one, much larger than one (Bi = 
10), and much lower than one (Bi = 0.1). For this specific system, Biot number Bi = 1 
correspond to initial height 1000 nm and evaporation rate k e = 0.001 m/s. 

When Biot number is much larger than one, evaporation is the dominant process. This 
is reflected in shorter total time of evaporation compared to other cases (see Figure [3] top) . 
Consequently, the solvent removed from the top layer cannot be balanced by mass diffusion 
within the film. Thus, a boundary layer lean in solvent and rich in other components is 
created close to the top (Figure [3] (right)). Enrichment of polymer and fullerene results in 
the initiation of phase separation. In this way, boundary layer becomes the region of the thin 
film when solution is unstable even to small fluctuations, which leads to phase separation. 
Phase separation subsequently propagates into the depth of the film. This is clearly seen in 
Figure [3] (right). Close to the top layer, a blocking layer rich in polymer is created. Because 
of low diffusion coefficient of polymer, the top layer blocks solvent evaporation from the top. 

When Biot number is much smaller than one, diffusion is the dominant process. Cor- 
respondingly, the total time of solvent removal is much longer. The system has more time 
to balance solvent lost from the top layer. Diffusion is not suppressed by fast solvent re- 
moval (as for high Biot numbers) and consequently no top boundary layer rich in polymer 
is created. There is no significant gradient of the solvent in the composition within the film. 
Consequently, phase separation in initiated homogeneously along the height. 

Biot number provides important insight into interplay between evaporation and diffusion. 
It can be used to link two types of morphology evolution: homogeneous across the film or 
initiated close to the top surface. For symmetric systems (N p = Nf) and low Biot number, 
solvent content within the film is homogeneous. For such a scenario, an assumption of 
homogeneous decrease of solvent volume fraction within the entire volume is valid. Such 
assumption was made in a recent work [57]. However, if one of the non-solvent components 
has larger molar volume, and consequently large Ni - which is the case in OSC - this 
assumption is invalid and evaporation must be included in the model explicitly. 

The evaporation rate affects not only the region where phase separation is initiated but 
also affects the average size of the separated phases. In Figure [4j we plot morphology 
evolution for three different evaporation rates characterized by Biot number 0.03, 0.3 and 3. 
The magnitude of the evaporation rate affects the total time of the process (as shown in 
Figure [3]). For lower evaporation rate, and Biot number, total process time is longer. Once 
phases separates, domains rich in each component have additional time to coarsen and create 
larger domain (see Figure [4] left). For higher evaporation rate, and Biot number, total time 
is shorter. Once phases separates they have very short time to coarsen (see Figure [4] right). 
This is due to the fact that solvent is removed from the system rapidly, and consequently 
the diffusion coefficient of the film decreases significantly leaving the morphology frozen. 

7.2. Influence of blend ratio 

The blend ratio between polymer and fullerene is a key system variable during fabrication 
of OSC, that is additionally relatively easy to manipulate. The blend ratio affects the type of 
morphology that develops during spinodal decomposition. Two basic classes of morphologies 
typically develop — percolated morphology and morphology with islands. Intuitively, the 
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Figure 3: (Top) Change in film height in time for three Biot numbers: Bi = 0.1, 1 and 10. (Bottom) Volume 
fraction profiles of polymer along film height for three different Biot numbers. 
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former is more suitable for OSC than the latter. This is because the latter structure has 
more islands not connected to relevant electrodes and cannot provide useful pathways for 
charges to reach the boundaries. Therefore, for OSC application, it is necessary to clearly 
identify process and system variables that lead to the percolated type of structure. 

We investigate two blend ratios 1:1 and 1:0.8. In Figure [5| we plot snapshots of the 
morphology evolution with time. Percolated morphology develops for blend ratio 1:1. How- 
ever, a small change in blend ratio of the analyzed system results in significantly different 
morphology type. This is clearly seen for blend ratio 1:0.8, where multiple stripes form spon- 
taneously and break up as the solvent evaporates. Such multiple layer formation was also 
observed in experiments and is reported [58]. We note that although a layered morphology 
has no application in OSCs, such morphology is desired in organic transistors [59j. 

7.3. Influence of degree of polymerization 

The degree of polymerization, N p) is another tunable variable that allows for controlling 
morphology evolution. In practice, for OSC fabrication, only polymer degree of polymeriza- 
tion can be controlled and it has been shown experimentally that efficiency is affected by 
its choice [60] . In Figure [6| we compare the morphology evolution for three different degree 
of polymerization: N p = 5, 20 and 100. All simulations have been performed for the same 
blend ratio 1:1 and Biot number Bi = 0.4. As we increase the degree of polymerization, 
the morphology changes significantly. For the symmetric case: N p = Nf = 5, a percolated 
morphology evolves; while for asymmetric cases multiple layers are created. Notice also that 
with increasing degree of polymerization, phase separation is initiated earlier. Moreover, 
when polymer degree of polymerization is larger than fullerene degree of polymerization, 
morphology type changes from percolated into multiple layered morphology. 

7.4- Influence of solvent type 

In solvent-based fabrication, the solvent creates the environment for morphology evolu- 
tion. Both components must be soluble in the common solvent to create an initially homo- 
geneous solution. Choice of the solvent has significant effect on morphology evolution and 
provides additional system variable for morphology control. This is manifested as differences 
in relative values of interaction parameters between all three components. 

In Figure [7| we show morphology evolution for three combinations of the interaction 
parameters (Xp/? Xps, Xfs)'- (1,0.3, 0.3), (1,0.3, 0.6) and (1,0.6, 0.3). In the first case, solvent 
is chosen such that interactions between solvent and two components are the same and 
much lower than interaction between polymer and fullerene. This means that polymer and 
fullerene have similar solubilities in solvent. In the second case, fullerene is less soluble in 
the solvent compared to polymer. Interaction parameter between fullerene and solvent is 
two times higher than between polymer and solvent. In the third case, we assume that 
polymer is less soluble in the solvent. In all three case, we assume the same evaporation rate 
(Bi = 0.3), blend ratio (1:1), and degree of polymerization (N p = 100, Nf = 5 and N s = 1). 

Changing solvent results in dramatically different morphology evolution [61]. Multiple 
layer formation that we observe for higher polymerization is broken in the third case, due to 
the different solvent used. 
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Figure 5: Morphology evolution for two different blend ratio between polymer and fullerene 1:1 and 1:0.8. 
Consecutive rows correspond to the morphology at height: initial, 0.6, 0.5, 0.4 and the final height. 
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Figure 6: Morphology evolution for three different degree of polymerization: N p = 5, N p = 20 and N p = 100. 
Consecutive rows correspond to the morphology at height: initial, h = 0.7, 0.6, 0.5, 0.42, 0.4, 0.35 and the 
final height. 
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Figure 7: Morphology evolution for three different solvent types characterized by different interaction pa- 
rameters between components. Consecutive rows correspond to the morphology at height: initial, h = 0.9, 
0.8, 0.7, 0.6, 0.5, 0.45 and the final height. 
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7 5. Effect of substrate 

The height of a typical OSC active layer is around 100 — 300 nm. For such thin ge- 
ometries, morphology close to the substrate can be significantly different than that of the 
bulk. We consider the effect of chemical patterning of substrate to selectively affect each 
component. We consider two cases. In the first case, the substrate is patterned to attract the 
polymer preferentially. In the second case, substrate is patterned with two chemistries: one 
preferentially attracting polymer and one preferentially attracting fullerene. In both cases, 
patterns of wavelength A s = 0.12 cover the substrate. In each case, we assume the chemical 
potential of patterned material ji p = \v = 0.01, and h p = h f = 0. In most cases, h is small 
compared to ji [62J. 

Figure [8] shows morphology evolution on three different patterned substrates: neutral, 
preferentially attracting polymer, and preferentially attracting polymer and fullerene. We 
run these tests for polymer of high degree of polymerization N p = 100 (last column in Fig- 
ure[6]right), which becomes the reference simulation and is repeated in Figure [8] (left). These 
results clearly demonstrate that substrate patterning provides additional degree of control over 
morphology. Multiple layers observed without patterning can be broken close to the sub- 
strate. The breakage depends on the frequency of the patterning and combination of the 
material types (Figure [8] middle and right). 

Patterning with alternating chemical preference allows for better control of domain size. 
The size along the substrate of the polymer-rich induced by such patterning is maintained 
during evaporation. This is not the case when patterning is purely of one preference. The 
domains created close to the substrate grow in size along the substrate. However, when 
substrate is patterned with one chemical preference, domains created close to the substrate 
penetrate deeper into film, compared to the other case. 

7. 6. Control of morphology evolution 

In previous subsections, we showed that by independently changing system and process 
variables we obtain various types of morphologies. In Figure [9| we summarize these types 
of morphologies. It is important to notice that significantly different morphologies can de- 
velop for various system and process variables. In general, there are three main types of 
morphologies: percolated morphology (b); morphology with multiple layers (c and d) and 
morphology with islands (e) . Substrate patterning gives an additional means to initiate and 
direct morphology close to the substrate. For example, adding substrate patterning leads to 
breaking the bottom layer, as shown in Figure [HJ It is interesting to note that controlling 
different variables may lead to the same type of morphology. For example, increasing degree 
of polymerization from N p = 5 to N p = 20 leads to multiple layer creation (Figure [6]). Similar 
effect is observed by changed solvent (Figure [7]). This emphasize the importance of further 
systematic studies of solvent-based fabrication. 

8. Conclusions 

Morphology is a key element affecting the performance of organic solar cells. The mor- 
phology evolution during solvent-based fabrication of organic solar cells is a complex, multi 
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Figure 8: Morphology evolution three systems where substrate is neutral to both components (left), substrate 
is patterned with material that attracts polymer (middle) and substrate is patterned with two materials, one 
attracting polymer, second attracting fullerene (right). Consecutive rows presents morphology evolution at 
h = 1, 0.7, 0.6, 0.5, 0.45 and final. 
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Figure 9: Morphology control by changing system and process variables of the solvent-based fabrication 
techniques. 
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physics process that is affected by a variety of material and process parameters. A vir- 
tual framework that can model 3-D morphology evolution during fabrication of OSC can 
relate these fabrication conditions with morphology. This will significantly augment organic 
photovoltaic research which has been predominantly based on experimental trial and error 
investigations. Such a framework will also allow high throughput analysis of the large phase 
space of processing parameters, thus yielding considerable insight into the process-structure- 
property relationships governing organic solar cell behavior that is currently in its infancy. 

In this work, we develop a phase field-based framework to study 3D nanomorphology 
evolution in the active layer of OSC during solvent-based fabrication process. In particular, 
we formulate physical and mathematical model that takes into account all important pro- 
cesses that occur during solvent-based fabrication of OSCs. We select phase field method 
to model the behavior multicomponent system with various driving forces for morphology 
evolution. We outline an efficient numerical implementation of the framework, to enable 
three dimensional analysis of the process. We showcase our framework by investigating the 
effect of various process and system variables, that lead to following observations: 

1. Mass Biot number expresses the interplay between solvent evaporation from the top 
surface and diffusion within the thin film. For high Biot number, evaporation is a 
dominant process which results in top boundary layer creation enriched in two main 
components. Consequently, phase separation initiates close to the top and propagates 
into the film. For low Biot number, in turn, diffusion is a dominant process, solvent 
removed from the top layer is diffused back from the resulting in homogeneous profile. 
Thanks to this uniform distribution phase separation initiates and evolve homoge- 
neously within film. 

2. The morphology evolution is affected not only by kinetics through evaporation and 
diffusion but also by thermodynamics. In particular, interaction parameters between 
components and degree of polymerization have a large effect on morphology evolution. 

3. The accessibility of possible configurations provided by the free energy landscape is 
controlled by system variable such as blend ratio. Small change of blend ratio lead to 
large variation in morphology evolution. 

4. Finally, surface induced phase separation provides another opportunity to locally affect 
the morphology, by creating additional "sinks" in the energy landscape. 

We are currently investigating extensions of the framework along three directions: nonho- 
mogeneous evaporation, fluid shear effects (based on [63]), and crystallization. 
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Figure A. 10: Size of time step profile and film height profiles for example ID, 2D, and 3D problems, respec- 
tively. 



Appendix A. Time step adaptivity 

We use Euler Backward scheme with the heuristic strategy to adjust size of time step 
This strategy is based on the number of Newton's iterations required to solve a nonlinear 
problem for given time step. If number of iteration is lower than 20 size of time step is 
increased by 25% in new time step, otherwise it is reduced to 25% of previous time step. 
When solution cannot be found in 50 iterations, such step is rejected and time step decreased 
by half. In this way time step is decreased when rare coarsening events occur and increased 



when morphology evolves slowly. In Figure |A.10[ we show example profile of time step size for 
various dimensionality considered. Size of time step is adjusted by few orders of magnitude. 
Efficient time stepping strategy allows to perform simulations that cover several time scales. 
In the figure, we show the time scale spanning over two to three orders of magnitude when 
height decreases by up to 80%. 
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